Iterative solution of a Dirac equation with inverse Hamiltonian method 
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We solve a singe-particle Dirac equation with Woods-Saxon potentials using an iterative method 
in the coordinate space representation. By maximizing the expectation value of the inverse of 
the Dirac Hamiltonian, this method avoids the variational collapse, in which an iterative solution 
dives into the Dirac sea. We demonstrate that this method works efficiently, reproducing the exact 
solutions of the Dirac equation. 
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The imaginary time method has been successfully em- 
ployed in non-relativistic self-consistent mean-field calcu- 
lations [l], [|| • The idea of this method is that a function 
e~ HT ' n \^)Q} converges to the ground state wave function 
of a Hamiltonian H as r — > oo for any trial wave func- 
tion \i[>o) as long as it is not an eigenstate of H. One can 
then iteratively seek for the eigenfunctions of H in the 
coordinate space starting from an arbitrary wave func- 
tion \ipo). This method is suitable particularly for self- 
consistent mean-field calculations in a three-dimensional 
coordinate space, with which an arbitrary shape of nuclei 
can be efficiently described 0. 

A naiive extension of this method to a relativistic equa- 
tion, however, meets a serious problem. That is, a Hamil- 
tonian H has both the Fermi and Dirac seas (see Fig. 
1) and an iterative solution inevitably ends up with a 
wave function in the Dirac sea even if the starting wave 
function |?/>o} is in the Fermi sea. This problem was 
recently pointed out in Ref. Q in the nuclear physics 
context, although the problem had been well known in 
quantum chemistry under the name of variational col- 
lapse [H-Qjl]. In Ref. [H, the problem was avoided by 
applying the imaginary time method to the Schrodinger- 
equivalent form of the Dirac equation. The same method 
was used also in Ref. [13| . Moreover, damped relaxation 
techniques were employed in Ref. [14j for an iterative 
solution of a Dirac equation. 

In this paper, we propose a yet another method for 
iterative solution of a Dirac equation, which can be di- 
rectly applied to the coordinate space representation of 
a relativistic wave function. The method is based on 
the idea of Hill and Krauthauser [8| to maximize the ex- 
pectation value of the inverse Hamiltonian. This method 
solves the Dirac equation as it is and thus both the upper 
and lower components of a wave function are automati- 
cally obtained. Also, in this method, states in the Dirac 
sea can be obtained within the same scheme by invert- 
ing the direction of "time" evolution, while the method 
in Refs. [1, [l3[ requires two different equations for the 
Fermi and Dirac seas. Thus, our method can be regarded 
complementary to the method in Refs. (3. Il3j. 

In our method, the ground state wave function of 
a Dirac Hamiltonian H is obtained by using the relation 
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FIG. 1: (Color online) (a) A schematic spectrum of a Dirac 
Hamiltonian H. V(r) and S(r) axe vector and scalar po- 
tentials, respectively. Eo is an arbitrary constant located be- 
tween the lowest energy state in the Fermi sea and the highest 
energy state in the Dirac sea. The hatched regions indicate 
the continuum spectra, (b) A spectrum of 1/{H — Eo) corre- 
sponding to Fig. 1(a). The filled circles indicate the bound 
states in the Fermi sea, while the open circles indicate the 
bound states in the Dirac sea. The continuum region is de- 
noted by the thick solid line. 



where Eq is a constant between the lowest eigenvalue in 
the Fermi sea and the highest eigenvalue in the Dirac 
sea. In order to illustrate how this method works, Fig. 1 
shows a spectrum of H (Fig. 1 (a)) and that oH/(H—Eq) 
(Fig. 1 (b)) H. In the spectrum of 1/(H - E ) shown 
in Fig. 1(b), the states in the Fermi sea appear in the 
positive region while those in the Dirac sea in the negative 
region. Therefore, as T — > oo in Eq. JT]), all the states 
in the Dirac sea are damped out. Moreover, the lowest 
energy state in the Fermi sea corresponds to the highest 
point in Fig. 1(b), and only this state survives in Eq. ([1]) 
as T — > oo. 

In practice, Eq. fTj is solved iteratively. That is, from 
the wave function at T = nAT, |\& n ), the wave function 
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is evolved with a small interval AT as 



with 
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The inverse of the Hamiltonian is in a familiar form 
in non-relativistic time-dependent Hartree-Fock (TDHF) 
calculations [TH, [01 (see also Ref. [HI for a similar tech- 
nique in relativistic TDHF calculations), and one can 
evaluate the inverse similarly to Refs. |15HL7j using the 
Gauss elimination method. We demonstrate it here for 
spherical potentials. The relativistic wave function then 
takes the form [ill [ill], 
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where I = 2j - I = j ± 1/2 for I = j =p 1/2, and 
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is the spin-angular wave function, Yi mi and Xm a being the 
spherical harmonics and the spin wave function, respec- 
tively. The Dirac equation with a spherical vector poten- 
tial V(r) and a scalar potential S(r) then reads [la. fl9|. 



he I 



U 



dr 



H-fr +f) 

W - 2mc 2 



f)= e (f ) ! "'> 



where k = =f(J+1/2) for j = l±l/2, U(r) = V(r) + S(r), 
and W(r) — V(r) — S(r). It is easy to show that G 
and F are proportional to r l+1 and r' +1 , respectively, 
around the origin, r ~ 0. Notice that |^ n+1 ) satisfies 
(i? - E )\f n+1 ) = (H - E + AT)|#") (see Eq. ©). 
This implies that G and T at T + AT satisfy 
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where G and F are defined as 
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We solve these equations by discretizing the radial coordi- 
nate with a spacing of Ar and imposing the box boundary 
condition, that is, the wave functions vanish at the edge 
of the box. Denoting the wave functions at r k = kAr 
as Gfc and F k , and using the three-point formula for the 
first derivative, i\)' h ~ (tpk+i — V'fc-i)/2Ar ) Eq. ([6]) reads 



4>k+\ + A k<t>k - 4>k-x 
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Here, A k is a 2x2 matrix defined by 
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We assume 



<t>k+l = a k<f>k+Pk, 



(10) 



(11) 



where a k is a 2x2 matrix and (3 k is a two-component 
vector. Substituting this to Eq. ([8]), one finds 



CKfc-l = (Afe + CCfe) 1 , 



(12) 
(13) 



These equations can be solved inwards from r^r-i with 
(Xn-i = = 0: which ensures that the wave func- 

tions vanish at i? max = fN- Once and /3 fe are so 
obtained, the wave functions <p k can be constructed with 
Eq. (TIT]) outwards from 4> k=Q = 0. 
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FIG. 2: The expectation value of the Hamiltonian H (Fig. 
2(a)) and that of the inverse of the Hamiltonian, 1/(H — Eq), 
(Fig. 2(b)) at each iteration for the neutron lpi/2 state in 
ls O. Eq is set to be the depth of the Woods-Saxon potential 
U(r) = V(r) + S(r), and the step AT is taken to be 10 MeV. 

Let us now numerically investigate the performance 
of the inverse Hamiltonian method. To this end, we 
use spherical Woods-Saxon potentials for U(r) and W(r) 
which correspond to 16 O. The parameters of the Woods- 
Saxon potentials are taken from Ref. [2(J. Figure 2 
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shows the convergence feature for the neutron lpi/2 state 
(« = 1), where the quantum numbers refer to the upper 
component of the wave function. The coordinate space 
is discretized up to i? max =15 fm with Ar = 0.1 fm. The 
energy shift Eq is taken to be the depth of the potential 
U(r), that is, E = -71.28 MeV, and the size of step AT 
is taken to be 10 MeV. For the initial wave functions, 
we take F(r) = and G(r) to be the non-relativistic 
lpi/2 wave function of a Woods-Saxon potential with 
V = -51 MeV, R = 1.27 x 16 1 / 3 fm, and a=0.67 fm 
21]. Fig. 2(a) indicates that the energy of the lpi/2 state 
is quickly converged as the number of iteration increases. 
The converged value is E = —18.974 MeV, that is almost 
identical to the exact value, E = -18.976 MeV. Fig. 2(b) 
also shows that the expectation value of the inverse of the 
Hamiltonian, 
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monotonically increases to the converged value, although 
it may not be the case for (H). 
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FIG. 3: (Color online) Comparison of the wave function for 
the neutron lpi/2 state in 16 O obtained with the inverse 
Hamiltonian method (the dashed lines) to the exact solutions 
(the solid lines). 

Figure 3 shows the wave functions for the lpi/2 state. 



The solid lines show the exact solution of the Dirac equa- 
tion obtained with the Runge-Kutta method, while the 
dashed lines show the wave function obtained with the in- 
verse Hamiltonian method. One can clearly see that the 
iterative method well reproduces the exact wave function. 

We have checked that the performance of the inverse 
Hamiltonian method is good also for other states with 
different quantum numbers as well, although the accu- 
racy seems to be somewhat improved if the energy shift 
Eq can be chosen as close as possible to the exact value. 
Evidently, the inverse Hamiltonian method provides an 
efficient method to solve a Dirac equation iteratively. 

Before we close this paper, we would like to point out 
that the states in the Dirac sea can be also obtained 
within the same scheme by changing AT by — AT in Eq. 
@. This is so because the highest energy state in the 
Dirac sea corresponds to the minimum of l/(i? — To) 
(see Fig. 1(b)). We have obtained the bound states in 
the Dirac sea in this way for 16 O, and have confirmed 
that the exact wave functions are well reproduced. 

In summary, we have proposed a new scheme to it- 
eratively solve a Dirac equation. The idea of the new 
scheme is to maximize the inverse of a shifted Hamilto- 
nian, 1/(H — Eq), and thus we call it the inverse Hamilto- 
nian method. By choosing Eq to be between the highest 
energy in the Dirac sea and the lowest energy in the Fermi 
sea, the expectation value of 1/(H — Eq) monotonically 
converges to the exact value, 1/(E — To). The upper 
and the lower components of a wave function can be si- 
multaneously obtained with this method in the coordi- 
nate space representation. We have applied this method 
to neutron states in 16 using spherical Woods-Saxon 
potentials. We have shown that the inverse Hamilto- 
nian method efficiently reproduces the exact energies and 
wave functions for the states both in the Fermi and Dirac 
seas. In this paper, we have assumed spherical symme- 
try for nuclear potentials. It will be an interesting future 
work to apply this method to self-consistent relativistic 
mean-field calculations in three-dimensional space. 
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